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Monte Carlo simulations are used to study lattice gases of particles with extended hard cores on a 
two dimensional square lattice. Exclusions of one and up to five nearest neighbors (NN) are consid- 
ered. These can be mapped onto hard squares of varying side length, A (in lattice units), tilted by 
some angle with respect to the original lattice. In agreement with earlier studies, the INN exclusion 
undergoes a continuous order- disorder transition in the Ising universality class. Surprisingly, we find 
that the lattice gas with exclusions of up to second nearest neighbors (2NN) also undergoes a contin- 
uous phase transition in the Ising universality class, while the Landau-Lifshitz theory predicts that 
this transition should be in the universality class of the XY model with cubic anisotropy. The lattice 
gas of 3NN exclusions is found to undergo a discontinuous order-disorder transition, in agreement 
with the earlier transfer matrix calculations and the Landau-Lifshitz theory. On the other hand, 
the gas of 4NN exclusions once again exhibits a continuous phase transition in the Ising universality 
class — contradicting the predictions of the Landau-Lifshitz theory. Finally, the lattice gas of 5NN 
exclusions is found to undergo a discontinuous phase transition. 

PACS numbers: 



I. INTRODUCTION 



Since the early days of statistical mechanics lattice 
gases serve as the foundation on which many models of 
complex physical systems are constructed. These range 
from simple fluids [lj to structural glasses and granular 
materials 0. The simplest lattice gas consists of non- 
interacting particles which are constrained to move on a 
lattice with a restriction that each lattice site is occupied 
by at most one particle. This model can be solved exactly, 
showing that the thermodynamics is trivial and no phase 
transition is present. A lattice gas in which particles in- 
teract with their nearest neighbors can be mapped onto 
the Ising model in external field and exhibits a first or- 
der phase transition terminating in a critical point. The 
difficulty and the scarcity of exact solutions has stimu- 
lated the development of various approximate theories 
aimed at treating more complex interaction potentials. 
Examples of these are the high and the low tempera- 
ture (or density) expansions generalization of Bethe's 
methods for Ising model (which became known as clus- 
ter variational methods) [J], as well as the approxima- 
tion schemes such as Rushebrook and Scoins yj. Since 
melting is dominated by strong short ranged repulsive 
forces, a lattice gas in which particles interact exclusively 
through an extended hard core — a particle on one site 
prevents the neighboring sites from being occupied - 
has attracted a particular attention Q . In these systems 
temperature plays no role since the interaction energy is 
infinite inside the exclusion region and vanishes outside. 
Surprisingly already with one nearest neighbor exclusion 
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(INN), an order-disorder transition appears. The tran- 
sition is purely entropic and is similar to the freezing of 
hard spheres, except that now the transition is of second 
order. 

In this paper we are interested in the phase transitions 
which occur in hardcore lattice gases as the size of the 
exclusion region is extended. We start with the INN 
and progressively increase the exclusion region up to 5 
nearest neighbors (5NN). To study the phase transitions, 
grand-canonical Monte Carlo simulations are performed 
with alternating attempts at inclusion/removal of par- 
ticles followed by some tentative diffusion moves. The 
grand-canonical ensemble was chosen because it allows 
for the density fluctuations which are particularly im- 
portant since the thermal fluctuations are absent in these 
systems. 

The specific question that we would like to address here 
is how does the order-disorder transition depends on the 
range of the exclusion region and whether the universal- 
ity class and the order of the phase transition can be 
predicted based purely on symmetry considerations. 

The paper is organized as follows: section [II] defines 
the relevant order parameters characterizing the ordered 
phases. Section IIIII presents the results of the grand- 
canonical Monte Carlo simulations, and section ITVl gives 
the final discussion and the conclusions. 



II. EXCLUSION MODELS, SUBLATTICES AND 
ORDER PARAMETERS 



A lattice gas is composed of particles whose positions 
are restricted to coincide with the vertices of a given lat- 
tice, and where each site can be occupied by, at most, 
one particle. When there is no further interaction be- 
tween the particles, this system presents a trivial ther- 
modynamics and no phase transition. Instead, here we 
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consider equilibrium lattice gases where the range of ex- 
clusion extends over the neighboring sites. Increasing 
the radius of the exclusion shell is equivalent to taking 
smaller lattice cell sizes as one goes to the continuum 
limit. These exclusion shells may be interpreted in terms 
of geometric hard core particles. All cases considered here 
are equivalent to parallel hard squares of size A, cither 
tilted or not. Since hard core potentials allow cither zero 
or infinite energies, temperature is not a relevant param- 
eter, the systems thus being athermal. As a consequence, 
the only relevant independent parameters are the volume 
(V) and the chemical potential (p), the latter determin- 
ing the average number of particles inside the volume. In 
the grand-canonical Monte Carlo simulation, three kinds 
of trial movements are allowed: displacement, insertion, 
and deletion of particles Q. The simulations are done 
in two stages. First various short runs are performed to 
locate the critical region. Then a single long run at the 
value of p close to the putative critical chemical potential 
p c is carried out. From the time series data of the long 
run, using histogram reweighting techniques d H fToj ] . 
the relevant thermodynamic obscrvablcs are accurately 
evaluated throughout the critical region. Another set of 
measurements at different values of p is also performed to 
check the reliability of the extrapolation. The diffusion 
trial movements are performed explicitly to help achieve 
higher densities and to to avoid entrapment of the system 
in low density metastable states. 

The lattice gases with hard core exclusion present a 
transition from a low density fluid-like phase, in which all 
sites have the same occupational probability, to a large 
density crystal-like phase in which the system is no longer 
translationally invariant. Thus, the lattice may be conve- 
niently divided into sublatticcs which in the high density 
phase are differently populated. Below we describe the 
exclusion shells considered in our simulations, with the 
corresponding sublattices. 

For the system sizes considered in our simulations, the 
density p and the compressibility k = L 2 ((p 2 ) — (p) 2 )/p 2 , 
where (...) is the ensemble average, are not very sensi- 
tive to the phase transition. To obtain a quantitative 
information we, therefore, build order parameters which 
become non zero as a system passes from a fluid to an 
ordered phase. When the transition is continuous, the 
order parameter q (to be defined below) obeys the usual 
finite size scaling relation 



q = L^^f[(p-p c )L 1 ^} 



(1) 



where / is a scaling function and p c is the critical value 
of the chemical potential in the thermodynamical limit, 
L — > oo. The staggered susceptibility x measures the 
fluctuations of q, 



X = L> ((, 2 ) - (qf) 



and has the scaling 

X = L^ v g[{p-p c )L l ' v ], 



(2) 



(3) 



where g is another scaling function. The exponents 
/?, 7 and v are related with the critical behavior of 
the order parameter, susceptibility and the correlation 
length, respectively, and obey the hyperscaling relation 
2(3/v + ^//i/ = d. In the simulations below, one way to 
obtain v is from the behavior of dlaq/dp, whose scaling, 
from cq. [1] follows 



d\nq 
dp 



L^ v h[{n - p c )L 1/h 



(4) 



where h(x) = f'(x)/f(x). An analogous equation is 
obeyed if we change g by x- Alternatively, v can also 
be obtained from the behavior of the shifted location of 
the maxima of x and k which approach the thermody- 
namic limit with L 1 ^ corrections. Once v is known, /3 
and 7 are obtained from eqs. Q] and El 



A. Nearest neighbor exclusion (INN) 

The simplest non trivial lattice gas consists of par- 
ticles which preclude nearest neighbor sites from being 
occupied. As shown in Fig. [T^, the lattice may be sub- 
divided into two sublattices, each site being surrounded 
by the sites of the other sublattice. This system, which 
can also interpreted either as 45° tilted hard-squares of 
linear size A = V2 or as hard disks of radius y/2/2, has 
been extensively studied and here we present some re- 
sults for the sake of both completeness and compari- 
son. Many different approaches have been used to de- 
scribe its properties on a square lattice: series expan- 
sions [!, El EH (3,03 (Cluster variational and transfer 
matrix methods [| Q GHEIEH El El HI EH El El , 
renormalization group 2 4~ 2 a| . Monte Carlo simula- 
tions!^ |2j| El EflEllI, E3, Bethe lattice @, 

MMm 

theory [3 



1371 . l3a ] , and more recently density functional 
Moreover, this model has also been consid- 
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ered because of its interesting mathematical 
and dynamical [H, 0, El Sfl SH El SI El 
properties. These studies indicate that the system un- 
dergoes a continuous phase transition belonging to the 
Ising universality class, at p c ~ 1.33, when the con- 
centration of particles is close to p = N/L 2 ~ 0.37. 
Note that the maximum density is 1/2. In higher dimcn- 
sionsEl [H EH, Hi El El, EH El and other geometries 
(see [40| and references therein) , the same kind of transi- 
tion is observed, also belonging to the Ising universality 
class. 

The transition is from the disordered fluid phase, where 
both sublattices are equally occupied, to a high density 
ordered phase with most of the particles confined to one 
sublattice (see Fig. []Ji). The order parameter for this 
system is 
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FIG. 1: Lattice division for the (a) INN and (b) 2NN exclu- 
sions: two and four sublattices are needed, respectively. The 
INN case is equivalent to tilted, non overlapping hard squares 
of length A = \/2, while in the 2NN case the squares are not 
tilted and have A = 2. The equivalent hard squares are shown 
in gray. 



where 



1,2 



(6) 



is the density in the ith sublattice Ci (the index j running 
over all sites in that sublattice), p M Ax = 1/2 is the maxi- 
mum equilibrium density. The variable aj is if the site 
j is empty and 1 if it is occupied. In the disordered low 
density regime pi = p%, and q\ = in the thermodynamic 
limit. At higher densities, one sublattice becomes prefer- 
entially occupied, resulting in a finite value of the order 
parameter. For finite systems the particles switch from 
one sublattice to the other, analogous to the switching 
of the magnetization sign in the Ising ferromagnet [59l |. 
The phase transition appears only in the thermodynamic 
limit. 



B. Nearest and next-nearest neighbors exclusion 
(2NN) 

Extending the range of interaction, we now consider 
a lattice gas such that each particle prevents both its 
nearest and next-nearest neighbors from being occupied. 
Thus, we define four different sublattices labeled from 
1 to 4 as shown in Fig. [Ux It should be noted that 
this system is equivalent to hard-squares with linear size 
equal to two lattice spaces (A = 2), as is demonstrated in 
Fig- El 3 ; if a particle occupies a position on one sublattice, 
say 1, the companion sites from sublattices 2 to 4 will be 
empty. Thus, the maximum possible density, correspond- 
ing to a full lattice occupation, is p M Ax = 1/4- Although 
this system has been considered many times in the last 

decades SJIsUp, HE IE iE [Ml III 11 HI, M EE 



7l|, sometimes as the zero temperature 
limit of a soft-core potential [U, [2(| [6(| [7l[ , there is still 
uncertainty of its universality class.. The order, location, 
and even the existence of the transition depends on the 
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TABLE I: Chemical potential and the density at the order- 
disorder transition along with the technique used and the pu- 
tative order of the transition for a lattice gas with nearest and 
next-nearest neighbors exclusion (2NN). The symbol "?" is 
used when some uncertainty is acknowledged by the authors. 
The several techniques used are: transfer matrix (TM), clus- 
ter variational method (CVM), low and high series expansion, 
generalizations of the Bethe method, interface method [74| . 
density functional theory (DFT) and Monte Carlo. Notice 
the wide range of values found for the critical chemical poten- 
tial. 



approximation method used to study the model. In ta- 
ble [J we show the density, the chemical potential, and 
the order of the transition obtained by several different 
methods. Little is known about its critical exponents, 
Ref. [2l|. 

In the low density regime, the system is fluid-like, with- 
out any long ranged order. As the density increases, there 
is a transition from the disordered phase to a columnar 
phase [1 i, HH, M, EE 111, where half-filled columns 
(rows) intercalate with empty ones: the system is or- 
dered along one direction but fluid along the other. Since 
there is no interaction between the particles (apart the 
exclusion), there is no alignment between particles in the 
neighboring columns (rows), that is, there is no solid 
phase (72)]. As a consequence, the close packed config- 
uration is not unique — each column (or row, but not 
both) may be shifted by one lattice space thus introduc- 
ing an O(L) contribution to the entropy at this maxi- 
mum density. Interestingly, no transition was found for 
the analogous d = 3 system of hard-cubes with A = 2 [|58[ 
on a cubic lattice. On the other hand, nearest and next 
nearest neighbors exclusion in d — 3 presents a weak first 
order transition 15811. Other lattices have also been con- 
sidered, see Refs. [40L [zl] and references therein. 

The order parameter needed to study this transition is 
a generalization of the INN Eq. O 
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In the fluid phase all the sublattices are equally occu- 
pied and q 2 = 0. The order parameter becomes non zero 
when the symmetry between the sublattices is broken. 
For the columnar phase two sublattices are preferentially 
occupied, while the remaining two have a much smaller 
density. For example, particles may be located on sub- 
lattices 1 and 2 (as in Fig. Q~p), while sublattices 3 and 4 
are almost empty etc. 



C. First, second and third neighbor exclusion 
(3NN) 

The 3NN case has its range of exclusion extended up to 
third nearest neighbors. As shown in Fig. [51 this can be 
interpreted both as hard-core pentamers [75| or slightly 
tilted squares with A = y/5. Previous studies of this 
system were based on series expansions la, ual . cluster 
variational and transfer matrix [f| [l8|, Wa, l77l. [78l . [79| , 
Bethe method etc.. It is now known that at high 
densities, the system undergoes a first order transition 
to a doubly degenerated ordered phase [8(|, the ground 
state configurations being related by chirality, see Figs. 
2b and 2c. This is reflected in the two possible ways of 
introducing sublattices to describe the symmetry break- 
ing: for each of the possible enantiomorph ground states, 
there is a corresponding labeling (called A and B) such 
that all particles belong to only one sublattice (see, for 
example, Fig. [5^ and b). 

For a given configuration, since we do not know a priori 
which of the two labcllings is more relevant, we measure 
the following quantity for each of them: 
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(8) 



and equivalcntly for the labeling B. Note that in the 
above sums the density of each sublattice appears four 
times, pi is given by eq. [5] for the chosen labeling. This 
quantity is zero for both labcllings in the fluid phase, 
where all sublattices are equally populated, while it be- 
comes non zero above the critical density (for one of the 
labcllings). The corresponding order parameter is the 
difference of the two measures: 

93 = — <I<7a - 9b |> • (9) 

Vmax 

Again, this order parameter is a natural extension of g 2 . 
At the close packed configuration, where the density is 
Pmax = 1/5, q A = 4p M Ax and q B = (or vice- versa): the 
first labeling is such that all particles sit on the same sub- 
lattice, while in the second labeling all sublattices have 
the same density. In both cases, however, (73 = 1. 



D. 4NN 

The subsequent situation considers exclusion of up to 
four nearest neighbor sites (4NN). As seen in Fig. [3Ji, 
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FIG. 2: Two possible labellings for the sublattices in the 3NN 
case — panels (a) and (b). For the same high density configu- 
ration in (a) all particles are on the same sublattice (5 in this 
example), while in (b) all sublattices are equally populated. 
For the labeling (b), panel (c) shows another configuration 
in which particles occupy only one sublattice. Note that the 
ground states shown in panels (b) and (c) are chiral, in the 
sense that the two are the reflections of one another about 
the left-to-right body diagonal. The exclusion problem can 
be formulated either in terms of the symmetric cross-shaped 
pentamers (shown in gray) or the tilted hard-squares of side 
length A = \fZ (d). 



this is equivalent to tilted hard squares of side length 
A = 2\/2. Nisbet and Farquhar [78|] studied this model 
through transfer matrices and concluded that the tran- 
sition is either continuous or very weak first order. The 
close packed configuration shown in Fig. is similar to 
the INN case, although less dense (the close packed den- 
sity is Pmax = 1/8), but it is not unique: analogously 
to the 2NN case, the columns (which in this case are 
tilted), are independent and may slide along the diago- 
nals. This does not change, however, the sublattice that 
the particles occupy. Notice also that this does not occur 
in the tilted diagonals of the INN and 3NN cases, the 
squares being always aligned. Thus, the order parame- 
ter should be insensitive to a diagonal displacement of a 
whole (tilted) column, and it is enough to subdivide the 
lattice in only two sublattices, in analogy with the INN 
case. Indeed, the transition to the ordered state may be 
analyzed with the same order parameter, cq. 
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FIG. 3: Sublattice division for the exclusion of up to (a) 4NN 
and (b) 5NN. These are equivalent to hard squares with A = 
2-\/2 and 3, respectively. Notice that in both cases the ordered 
phase is columnar (as in the 2NN case) , but in the former, the 
columns are along the diagonals. 



E. 5NN (A = 3) 
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FIG. 4: Density p as a function of the chemical potential fi 
for the INN exclusion for several lattice sizes. The inflection 
close to the transition point (/x c ~ 1.33) is very small, being 
noticeable only for larger system sizes. Inset: compressibility 
k as a function of the chemical potential. 



Finally, we considered a lattice gas composed of hard- 
squares with linear size equal to three lattice spaces, A = 
3. In terms of exclusion shells, this is also equivalent 
to excluding up to the fifth nearest neighbor (5NN), as 
shown in fig. [3Jd. The close packed configuration, like in 
the 2NN (A = 2) case, is not uniquely defined and has a 
density p MAX = 1/9. In d = 3, this system undergoes a 
weak first order phase transition as density increases [5c| , 
although there seems to be no information concerning 
the d = 2 case. The corresponding sublattice division is 
analogous to the 2NN case, only that now nine sublattices 
are necessary, the order parameter being 

95 = 7T- — (|Pi ~Ps\ + \pi ~Pq\ + \ps ~Pa\ 

+\P2 - Pel + \P2 - P7\ + \pe - py\ (10) 

+ |P3 - Pi\ + \P3 - Ps\ + \P4 - Ps\) (11) 

Eq. (fTT|) uses the labeling of sublattices depicted in 
Fig. (2,). 
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FIG. 5: Order parameter qi, Eq.[5] as a function of the chem- 
ical potential /i, for the INN exclusion for several lattice sizes, 
along with the corresponding data collapse (inset). The crit- 
ical exponents used are those of the two dimensional Ising 
model, v = 1 and /3 = 1/8. 



III. RESULTS 
A. INN 

In the case of nearest neighbors exclusion, the change 
in density at the transition is quite subtle, the inflection 
point being only noticeable for the largest sizes simu- 
lated, as seen in Fig. [H The transition is more clearly 
seen in the plot of the INN order parameter, Eq. [3 as 
shown in Fig. for increasing lattice sizes: q\ changes 
from a small value (fluid phase) to a value close to unity 
(checkerboard pattern). The behavior of qi is character- 
istic of a second order phase transition and indeed, all 
curves collapse onto a universal curve f\ , as shown in the 



inset of Fig. [5l whose scaling obeys eq. [TJ In accordance 
with the previous studies [H, [46|, this transition occurs 
at /i c ~ 1.33, corresponding to p c ~ 0.37, and belongs to 
the universality class of the two dimensional Ising model, 
with exponents v = 1 and (3 = 1/8. 

Another quantity of interest, shown in Fig. is the 
staggered susceptibility \i = ^ 2 ((<Zi) — (<7i) 2 ) which mea- 
sures the order parameter fluctuations. As the system 
size increases, the susceptibility peak becomes higher and 
narrower, shifting to larger values of the chemical poten- 
tial. All curves can be collapsed onto a universal curve 
whose scaling obeys eq. EH as shown in the inset of Fig. [5] 
Again, the exponents are those of the two dimensional 
Ising model, 7 = 7/4 and v — 1. Much less clear are the 
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FIG. 6: Order parameter fluctuations xij Eq. O as a function 
of the chemical potential /j, for the INN model using several 
lattice sizes. In the inset, the corresponding data collapse is 
shown using the critical exponents of the bidimensional Ising 
model, 7 = 7/4 and v = 1. 



FIG. 7: Density as a function of chemical potential for several 
lattice sizes for the gas of 2NN exclusion — region where 
the inflection point appears (noticeable only for the larger 
system sizes). Inset: compressibility n in the region where 
the transition is found. Only for the largest lattice simulated 
the maximum is prominent. 



fluctuations in density, that is, the compressibility k , as 
shown in the inset of Fig. [4j Although the peak grows 
with L, the sizes considered are still too small to try to 
obtain useful information from the data collapse. Notice 
that in this case, the maximum is expected to increase 
with the logarithm of L (a = 0). 

B. 2NN 

Analogously to the previous case, when the range of ex- 
clusion increases to include the next nearest neighbors, it 
is very hard to see the transition by only looking at the 
behavior of the density as a function of the chemical po- 
tential. In the critical region only for larger system sizes 
the inflection point becomes noticeable, as is shown in 
Fig. [JJ and the change in density is roughly one order of 
magnitude smaller as compared with the INN case. On 
the other hand, the order parameter q2 , Eq. shows a 
fast increase between \x = 4 and 5, Fig. [51 as the system 
enters the ordered (columnar) phase from the disordered 
(fluid) one. The the existence, location, order, and ex- 
ponents of this transition have remained elusive, as can 
be seen from a very broad range of approximate values 
presented in tableland the lack of any estimates for the 
critical exponents. We find that neither the density nor 
the order parameter q^ present any signs of discontinuity. 
Within the range of lattice sizes considered no hystere- 
sis, discontinuities, nor any evidences of the coexistence 
in the density and or the order parameter histograms has 
been observed. We conclude that the transition is con- 
tinuous, in agreement with most of the previous studies. 

Indeed, the q^ data can be collapsed onto a universal 
curve (as in Eq. Q}, with [i c ~ 4.574, f3/v = 0.125 and 
v = 0.94. Remarkably, j3Jv has the same value of the 
(exact) Ising model, 1/8 = 0.125. Although v is a lit- 




3.5 4 4.5 5 5.5 6 



[' 

FIG. 8: Order parameter for 2NN (eq.[7J as a function of the 
chemical potential for different lattice sizes. Small values of 
q2 signal that the system is disordered. For qi closer to unity, 
the system is in a columnar phase (see text). Inset: data 
collapse onto a universal curve for the larger system sizes L = 
128, 160, 256, 320 and 360. A very good collapse is obtained 
for jj, c = 4.574, v — 0.94 and j3/v = 0.125, in close agreement 
with the exact Ising values. The collapse remains good if we 
use the Ising value v = 1. 



tie bit smaller than the Ising value, a still good collapse 
is obtained with v = 1 and, as we will sec below, is in 
agreement with other measures of this exponent. The 
density, measured for several systems sizes at fi Cl or at 
the peak of the susceptibility extrapolate when L — > oo 
to p c = 0.233. The value of fj, c is also compatible with 
the crossing point of the Binder cumulant for the order 
parameter and with the extrapolated position of the sus- 
ceptibility, fig. M Mc - 4.578. 
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FIG. 9: Several estimates of the exponent v (2NN): the po- 
sition of the maxima of eq. [4] (top curve) and the suscepti- 
bility \2 (bottom), that shift as L~ 1//l/ , as a function of L. 
The solid lines are power-law fits neglecting the smaller sizes, 
from which we get v ~ 1.042 and v ~ 0.968, top and bottom, 
respectively. These values are very close to the Ising v = 1, 
as can be seen by the fits fixing the exponent to this value 
(dashed lines). Also from the fit of the location of the maxi- 
mum of the susceptibility we have an independent estimate of 
the transition: fi c ~ 4.578. Inset: the height of the maximum 
of eq. [4] increasing as L 1 ^ , as a function of L. Again the 
dashed line is a fit with v = 1, while the solid line is the best 
fit with v ~ 0.937. 



We also measured the order parameter and the density 
fluctuations — that is, the staggered susceptibility X2, 
eq. [2J and the compressibility «2, respectively, as shown 
in Figs. [TO] and [7] (inset) as a function of the chemical 
potential. In both cases, as the system size increases the 
curves get less broad and their height becomes larger. 
However, the maximum of K2 only becomes noticeable for 
the larger simulated lattices and no reliable information 
can be obtained from its scaling, within the range of L 
considered here. The shift of the position of the peak 
of X2 from (j, c behaves as L^ x l v (in the same way as 
the position of the maximum in eq. [4]), from where v 
can be estimated, as shown in fig. [5J In addition, the 
height of the peak of eq. d] increases with L 1 ^ . The 
exponent v evaluated from both measurements is very 
close to the Ising value, see caption of fig. [9] for details. 
We conjecture, therefore, that the exact value is indeed 
v = 1. The height of the peak of \2 increases as L 1 ^, and 
the fit yields j/v = 1.755 (bottom inset of fig. [TO)) , giving 
an excellent data collapse as can be seen in the top inset 
of fig. [TOJ Again, this value is very close to the known 
exact value for the Ising model, 7 = 7/4 = 1.75. The 
fact that the lattice gas of 2NN exclusion is in the same 
universality class as the Ising model is quite remarkable 
considering that the symmetries of the ground state and 
of the order parameter are so very different from those of 
INN model. 




FIG. 10: The staggered susceptibility X2 as a function of \i 
for several lattice sizes in the 2NN case. Inset: data collapse 
(top) onto a universal curve with exponent 7 obtained from 
fitting the height of the maximum of \2 as a function of L 
(bottom), 7/1/ = 1.755, along with v = 1. 



C. 3NN 

The lattice gas with exclusion of up to three nearest 
neighbors (3NN) is important since it presents, differently 
from the previous cases, a first order transition, showing 
that the nature of the transition depends on the range of 
exclusion, as was observed previously. Recently, Eisen- 
berg and Baram (7fjl | precisely located the transition, us- 
ing matrix methods, obtaining [i c ~ 3.6762. In Fig. [TT1 
and its inset, both the density and the order parameter, 
respectively, show a rapid rise as the chemical potential 
increases, passing from a fluid, low density phase, to an 
ordered, high density phase. Notice also the good agree- 
ment between the matrix prediction for the transition 
and our Monte Carlo data. The first order character of 
this transition is clearly seen in the behavior of the com- 
pressibility [7fj| and the fluctuations of the order param- 
eter, shown in Figs. [TOl and [TOl All the relevant scaling is 
with the volume of the system L 2 (y = 1/d [ll[), as can 
be seen in the inset of both figures. Excellent data col- 
lapses were obtained with fj, c = 3.6746, very close to the 
value calculated using the transfer matrix [79( . This value 
is also very close to the ones obtained by extrapolating 
the positions of the compressibility and the susceptibil- 
ity peaks as a function of the system size, and with the 
crossing points of the density and the order parameter 
for the two largest system sizes. 



D. 4NN 

When excluding up to four neighbors, the close packed 
state resembles both the INN and 2NN case. On one 
hand, particles are located on one of the two sublattices, 
arranged along the diagonals of the system (although 
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FIG. 11: Density as a function of chemical potential for the 
3NN exclusion model. At low densities the system behaves 
as a fluid and undergoes a first order phase transition to an 
ordered phase as the density increases. The vertical dashed 
line at /j, — 3.6762 is the infinite size extrapolation obtained 
with matrix methods [z^] and agrees well with the crossing 
point of the curves. Inset: the order parameter as a function 
of the chemical potential. 
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FIG. 12: Compressibility k as a function of the chemical po- 
tential for the 3NN exclusion for several lattice sizes. Inset: 
after rescaling the height and the width of the curves by L 2 
and L~ 2 , respectively, an excellent collapse is observed for the 
two largest systems with /i c = 3.6746. 



skipping one every two sites); on the other hand, these 
diagonals arc independent "columns" that may slide, orig- 
inating a 2 L -degenerated ground state with density 1/8. 
Surprisingly, within the size range available, the transi- 
tion appears to be continuous and in the Ising universal- 
ity class, contradicting the Landau-Lifshitz theory which 
predicts that this transition should be of first order [64 1 . 

As in the 2NN case, several estimates of v can be ob- 
tained. From the positions of the maxima of eq. 3] and 
Xi we obtain, fig. Q31 v ~ 0.923 and v ~ 0.999, respec- 
tively. The last one is more reliable since the peaks of 
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FIG. 13: Staggered susceptibility \3 as a function of the chem- 
ical potential for the 3NN exclusion case for several lattice 
sizes. Inset: after rescaling the height and the width of the 
curves by L 2 and L~ 2 , respectively, an excellent collapse is 
observed with [i c = 3.6746. 



Xi are less broad, thus allowing a more precise location 
of the maxima. A further measure of v is obtained from 
the height of eq. Hand gives v ~ 1.029. As in the 2NN 
case, these values oscillate around the Ising v = 1 and, 
again, we conjecture that this is the exact value. Indeed, 
this exponent also fits quite well the data. As before, 
the transition is not clear from the behavior of the den- 
sity, while a properly defined order parameter is able to 
capture it well, fig. 1151 The transition from the ordered, 
tilted columnar phase, to the fluid phase is found to occur 
at fji c ~ 4.705, corresponding to p c ~ 0.110. 

The order parameter and its collapse are shown in 
fig. [T5J The collapse is not very sensitive to the chosen 
value of (3, and any value in the broad range between 0.12 
and 0.14 will do well. We thus present the collapse with 
P/v = 0.125, the Ising value. Although this is rather ar- 
bitrary, the closeness of both v and 7 (see below) to the 
Ising values inspire us to believe that indeed this system 
is in the Ising universality class. 

The fluctuations of the order parameter are shown in 
fig. HU The dependence of the susceptibility peak on 
the system size is shown in the inset (below), yielding 
7/V = 1.681. With these values for the exponents, all 
curves can be collapse onto a universal one. An equally 
good collapse is obtained if instead of this value of 7 
we use the Ising one, 1.75. We conclude that the phase 
transition in the lattice gas with 4NN exclusion is in the 
Ising universality class. 



E. 5NN (A = 3) 

The lattice gas with exclusion of up to five nearest 
neighbors is equivalent to hard square particles with side 
length A = 3. In this case, as for the 2NN and 4NN, the 
system does not have a uniquely defined close packed con- 
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FIG. 14: Several estimates of the exponent v (4NN): the posi- 
tion of the maxima of eq.[4]and the susceptibility xi (bottom), 
that shift as L _1// ", as a function of L. The solid lines are 
power-law fits neglecting the smaller sizes, from which we get 
v ~ 0.923 and v ~ 0.999, top and bottom curves, respectively. 
These values are indeed very close to the Ising one, v = 1 (no- 
tice also that the peaks of eq.[4]are broader than those of xi, 
and their positions are less reliable). Also from the fit of the 
location of the maximum of the susceptibility we have an in- 
dependent estimate of the transition point: /j, c ~ 4.705. Inset: 
the height of the maximum of eq. [4] increasing as L 1 '" , as a 
function of L. The solid line is a fit with v ~ 1.029. 
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FIG. 15: Order parameter for the 4NN exclusion (eq. [5} as 
a function of the chemical potential for different lattice sizes. 
For large values of qi , the system is in a columnar phase with 
columns aligned along the diagonals. Inset: collapse of data 
onto a universal curve with v = 1, j3 = 0.125 and /i c = 4.705. 
Equally good collapse is obtained with a broad range of /3. 

figuration which can bring problems when approximated 
(analytical) methods such as series expansions and clus- 
ter variational methods are applied. 

When the density is measured as a function of the 
chemical potential, fig. [TTJ there is a fast increase be- 
tween [i = 5.5 and 5.6, and the curves become steeper as 
the system size gets larger, signaling the ordering of the 
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FIG. 16: The staggered susceptibility xi as a function of \i 
for several lattice sizes for 4NN exclusion. Insets: (top) data 
collapse onto a universal curve with exponent 7 obtained from 
fitting the height of the maximum of xi as a function of L 
shown in the bottom inset; j/v = 1.681 and v = 1. 



system. However, the increase in the density is rather 
small, while a much more pronounced change is observed 
for the order parameter, eq. [TTJ The fluctuations of the 
order parameter, shown in Fig. 1181 can be well collapsed 
with a L~ 2 scaling (the volume of the system), at least 
for the larger sizes considered. This is a signature of a 
first order transition. The location of the maximum de- 
pends on L and obeys the L~ 2 shift with fi c ~ 5.554. 
There seem to be, however, strong finite size effects for 
smaller system sizes, as observed when trying to collapse 
all the curves — only for the larger sizes considered the 
collapse is sufficiently good. For the density fluctuations 
(not shown), even stronger effects are seen. This is prob- 
ably due to the fact that the increase in the compress- 
ibility with L is very small and the non singular part of 
the free energy can not be discarded when compared to 
the singular one. This seems to be true for all the cases 
considered in this paper (being least pronounced for 3NN 
exclusion). 



IV. CONCLUSIONS 

We used the grand-canonical Monte Carlo simulations 
to study the order-disorder transition in lattice gases 
composed of particles with extended hard core on a two 
dimensional square lattice. In agreement with earlier 
studies the INN exclusion — corresponding to hard 
squares of length y2 tilted by 45° — was found to un- 
dergo a continuous order-disorder transition in the Ising 
universality class. The lattice gas with exclusions of up 
to second nearest neighbors (2NN) was also found to un- 
dergo a continuous phase transition in the Ising univer- 
sality class, while the Landau-Lifshitz (LL) theory [64| 
predicts that this transition should be in the universality 
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FIG. 17: Density as a function of the chemical potential for 
5NN exclusion (A = 3) for different lattice sizes. The steep- 
ness increases as the system develops a discontinuity. Inset: 
order parameter q$ as a function of /i for the same system 
sizes. 
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FIG. 18: Staggered susceptibility \5 as a function of the chem- 
ical potential for the 5NN exclusion (A = 3) for several lattice 
sizes. Inset: after rescaling the height and the width of the 
curves by L 2 and L~ 2 , respectively, a good collapse is ob- 
served for the larger simulated sizes with fi c — 5.554 



class of the XY model with cubic anisotropy. The lattice 
gas with exclusions of up to 3NN was found to undergo 
a discontinuous order-disorder transition, i n ag reement 
with the earlier transfer matrix calculations 17911 and the 



LL theory [64[. On the other hand, the gas of 4NN exclu- 
sions once again exhibited a continuous phase transition 
in the Ising universality class — a conclusion which con- 
tradicts the LL prediction of a first order phase transi- 
tion ■ Finally, the lattice gas with exclusions of up to 
5NN was found to undergo a discontinuous phase transi- 
tion, in agreement with the LL theory. 

The failure of the symmetry based Landau-Lifshitz 
theory to correctly predict the universality class and the 
order of the phase transition for some two dimensional 



lattice gases is not very surprising. It has been known 
for a long time that while the Landau-Lifshitz theory 
works well in three dimensions, its foundations are much 
shakier in 2d. For example, one of its requirements for 
the existence of a continuous phase transition is the ab- 
sence of third order invariants of the irreducible repre- 
sentation of the space group of the high symmetry phase 
to which the order parameter belongs. In 2d this condi- 
tion is known to break down — both three and four state 
Potts models undergo a continuous phase transition, even 
though their Landau- Ginzburg- Wilson (LGW) Hamilto- 
nians possesses third order invariants [821 . |83| . In view 
of this problem, Domany et al. have avoided the use of 
this "Landau rule" in their classification of the continuous 
order-disorder transitions of 2d lattice gases [64[ . Never- 
theless, the fact that one of the Landau rules fails, can 
already serve as a strong warning against a complete re- 
liance on the symmetry arguments for two dimensional 
systems. The fundamental difficulty is that 2d is the 
lower critical dimension for lattice gases. The fluctua- 
tion are very important and can easily modify the or- 
der of the phase transitions from the one predicted by 
mean-field. Furthermore, in 2d one can no longer rely on 
the renormalization group theory (RG) to justify a sim- 
ple low order polynomial form of the LGW Hamiltonian. 
Unlike in 3d, where it is possible to show that the high 
order invariants in the LGW Hamiltonian renormalize to 
zero under the action of RG, in 2d all the operators are 
equally relevant — thus invalidating a simple low order 
polynomial form of the LGW Hamiltonian. In the case 
of the isotropic XY model the situation is particularly 
dire since the symmetry based construction of the LGW 
Hamiltonian fails to consider the topological excitations 
which lead to the Kosterlitz-Thouless phase transition in 
this system [HI, [8^, [Hj]. It is, therefore, unlikely that 
symmetry considerations alone will be sufficient to pre- 
dict the universality class and the transition order for all 
2d lattice gases. 

Nevertheless, while the second order transition for the 
gas of 4NN exclusions is completely incompatible with 
even the modified version of the LL theory (in which the 
third order invariants are allowed), the case of 2NN is not 
nearly as bad. The LGW Hamiltonian for the lattice gas 
with 2NN exclusion [H is 



H = 



d 2 x 



{ (V^) 2 + (V^ 2 ) 2 +r(^ + ,/; 2 2 ) + 



(12) 



where ^ with i = 1,2 are the two components of the 
order parameter belonging to the two dimensional irre- 
ducible representation of the space group P4mm. 

In the absence of cubic anisotropy (last term of 
cq. (|12p ) this Hamiltonian will not have a phase transi- 
tion since the massless Goldstone excitations will destroy 
the long-range order at any finite temperature. Further- 
more, cq. (fT2")l does not account for the topological defects 
(vortices) which are responsible for the phase transition 
in the real XY model 0, [85|, HH and the appearance of 
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a pseudo-long-range order. Thus, in the passage from a 
microscopic Hamiltonian to the coarse grained LGW de- 
scription, important physics has been left behind [87ll88j|. 
Clearly under such conditions no conclusions based on 
this LGW Hamiltonian can be trusted. Nevertheless, if 
we take the form of eq. (fT2|) more seriously than per- 
haps it deserves, the Ising criticality is not incompatible 
with its structure. Based on the general properties of 
RG flows, we expect that eq. (fT2|) will have an Ising 
fixed point at which the isotropic fourth order invariant 
will vanish. At this fixed point, the LGW Hamiltonian 
will decouple into two Ising-like Hamiltonians for each 
component of the order parameter. The symmetry ar- 
guments, however, cannot tell us if the Ising fixed point 
is stable, or even if it is accessible to the renormaliza- 
tion group flow. Nevertheless, unlike the second order 
phase transition for 4NN exclusions, the Ising criticality 
for 2NN exclusions does not, in principle, invalidate the 
predictions of the LL theory. This not withstanding, it 
is quite surprising that in a multidimensional parame- 
ter space with an infinite number of fixed points (fixed 
lines) known to exist for the real XY model with cubic 
anisotropy 89], the 2NN model just happens to be in the 
region dominated by the Ising fixed point. 

For lattice gas with 4NN exclusion the simulations find 
a continuous order-disorder transition in the Ising univer- 
sality class. We note, however, that within a simulation it 
is impossible to be absolutely certain that the transition 
is not a very weakly first order as happens, for exam- 
ple, for the five states Potts model. There is, however, 
a big difference between q states Potts models and hard 
core lattice gases studied in this paper. For Potts model, 
the mean-field theory [83[ predicts that all transitions 
for q > 3 are of first order. The fluctuations, however, 
modify the order of q = 3 and 4 to second, and turn 
q = 5 very weakly first order [HJ . Importance of fluctu- 
ations decreases for larger values of q and the mean-field 
becomes progressively more accurate — the phase tran- 
sitions for q > 4 are all of first order. For lattice gases, 
on the other hand, mean-field theory (LL) predicts the 
transition to be of second order up to 2NN exclusions and 
first order afterwards. Indeed, we find continuous phase 
transitions for INN and 2NN exclusions and a first order 
transition for 3NN. Unlike the case of Potts model, how- 
ever, the agreement with the mean-field theory is broken 
precisely at 4NN and re-established again for 5NN. There 
are no precursor fluctuation induced second order phase 
transitions (such as q — 3 and 4 for Potts model) which 
would account for 4NN being weakly first order instead of 
second order. Furthermore, the critical exponents place 
4NN squarely into the Ising universality class, so that it 
would seem like a lot of coincidence for the Ising critical- 



ity to be only a crossover on the way to an underlying 
first order phase transition. We conjecture that, in fact, 
4NN is not the only continuous transition for higher order 
exclusions. This conjecture is based on the observation 
that as the exclusion region grows, the system approaches 
a continuum limit — lattice becomes irrelevant. For con- 
tinuum 2d systems melting, however, is driven by the 
topological defects (dislocations) and is of second (actu- 
ally infinite) order, belonging to the Kosterlitz-Thouless 
universality class [H [H, Hfl IS HH, Clearly, this 

can not be reconciled with the LL prediction of a first or- 
der transition for all exclusions above 2NN. Thus, there 
must be other lattice gases with exclusions larger than 
5NN which will also exhibit a continuous phase transi- 
tion. The only way to fit this into LL theory is to relax 
the so called "Lifshitz rule" which forbids the represen- 
tations for which the antisymmetric part of their direct 
product contains a vector representation. This, however, 
would completely destroy the LL theory since the num- 
ber of possible representations which can then be used 
to construct an invariant LGW Hamiltonian will become 
infinite. 

It is very difficult to understand the mechanism which 
changes a phase transition from first order, predicted by 
a mean-field theory, to second order when the fluctua- 
tions are taken into account. In two dimensions, how- 
ever, this seems to be a very common occurrence as is 
demonstrated by the 3 and 4 state Potts models which 
undergo continuous transitions in spite of third order in- 
variants present in their mean-field free energies. In this 
work we find that another mechanism is also possible 
— the fluctuations can invalidate the Lifshitz condition, 
transforming a first order phase transition into a fluctua- 
tion induced second order one. In view of our conjecture 
that INN, 2NN and 4NN are not the only lattice gases 
which undergo a continuous phase transition, it will be 
very curious to see what other higher order exclusions 
will also exhibit a continuous transition and what dis- 
tinguishes these systems from the ones which undergo a 
discontinuous order-disorder phase transition. 
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